//
//---------------------  1D  ---------------------
//
//
// Upwind is not defined for 1D
//

//
//---------------------  3D  ---------------------
//
//       Staggered Mesh for w-vel and v-vel
//
//   0       1       2       3       4       5   
//
//5      >       >       >       >       > 
//       |       |       |       |       |              
//   ^---+---^---+---^---+---^---+---^---+---^  4       Mesh for scalar fields
//       |       |       |       |       | 
//4      >   o   >   o   >   o   >   o   >               5  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |               4  x o | o | o | o x
//   ^---+---^---+---^---+---^---+---^---+---^  3           +---+---+---+---+
//       |       |       |       |       |               3  x o | o | o | o x
//3      >   o   >   o   >   o   >   o   >                  +---+---+---+---+
//       |       |       |       |       |               2  x o | o | o | o x
//   ^---+---^---+---^---+---^---+---^---+---^  2           +---+---+---+---+
//       |       |       |       |       |               1  x o | o | o | o x
//2      >   o   >   o   >   o   >   o   >               0  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |                  
//   ^---+---^---+---^---+---^---+---^---+---^  1           0 1   2   3   4 5 
//       |       |       |       |       |                 
//1      >   o   >   o   >   o   >   o   >                 o central node
//       |       |       |       |       |                 x boundary node
//   ^---+---^---+---^---+---^---+---^---+---^  0          > w velocity 
//       |       |       |       |       |                 ^ v velocity
//0      >       >       >       >       >    
//       0       1       2       3       4             
//
//                                                  
//          Volumes for w-velocity
//
//   0       1       2       3       4       5
//
//5      >       >       >       >       > 
//       :       :       :       :       :              
//   ^...+---^---+---^---+---^---+---^---+...^  4       Mesh for scalar fields
//       |   |   :   |   :   |   :   |   | 
//4      >   o   >   o   >   o   >   o   >               5  x-x-+-x-+-x-+-x-x
//       |   |   :   |   :   |   :   |   |               4  x o | o | o | o x
//   ^...+---^---+---^---+---^---+---^---+...^  3           +---+---+---+---+
//       |   |   :   |   :   |   :   |   |               3  x o | o | o | o x
//3      >   o   >   o   >   o   >   o   >                  +---+---+---+---+
//       |   |   :   |   :   |   :   |   |               2  x o | o | o | o x
//   ^...+---^---+---^---+---^---+---^---+...^  2           +---+---+---+---+
//       |   |   :   |   :   |   :   |   |               1  x o | o | o | o x
//2      >   o   >   o   >   o   >   o   >               0  x-x-+-x-+-x-+-x-x
//       |   |   :   |   :   |   :   |   |                  
//   ^...+---^---+---^---+---^---+---^---+...^  1           0 1   2   3   4 5 
//       |   |   :   |   :   |   :   |   |                 
//1      >   o   >   o   >   o   >   o   >                 o central node
//       |   |   :   |   :   |   :   |   |                 x boundary node
//   ^...+---^---+---^---+---^---+---^---+...^  0          > w velocity 
//       :       :       :       :       :                 ^ v velocity
//0      >       >       >       >       >    
//
//       0       1       2       3       4              
//
//                       
//                  |           |           |           |
//                --^-----------^-----------^-----------^-- 
//                  |     :     |     :     |     :     |               
//                  |     :     | (i,j+1,k) |     :     |
//                  o     >     o    w_N    o     >     o   
//                  |     :     |     :     |     :     |
//                  |     :     |     :     |     :     |
//                --^---------- 3 -- v_n -- 4 ----------^--  4 = v(i+1, j  ,k)
//                  |     :     |     :     |     :     |    3 = v(i  , j  ,k)
//                  |     :     |     :     |     :     |    2 = v(i+1, j-1,k)
//                  o    w_B   w_b   w_P   w_f   w_F    o    1 = v(i  , j-1,k)
//                  | (i,j,k-1) |  (i,j,k)  | (i,j,k+1) |
//                  |     :     |     :     |     :     |
//                --^---------- 1 -- v_s -- 2 ----------^-- 
//                  |     :     |     :     |     :     |               
//                  |     :     |     :     |     :     |
//                  o     >     o    w_S    o     >     o   
//                  |     :     | (i,j-1,k) |     :     |
//                  |     :     |     :     |     :     |
//                --^-----------^-----------^-----------^--
//                  |           |           |           | 
//                   
//   w_b = ( w(i,j,k-1) + w(i,j,k) ) / 2     
//   w_f = ( w(i,j,k+1) + w(i,j,k) ) / 2
//   v_n = ( v(i,j,k) + v(i+1,j,k) ) / 2   
//              3           4
//   v_s = ( v(i,j-1,k) + v(i+1,j-1,k) ) / 2
//              1           2

//   w_e = ( u(i,j,k) + u(i,j+1,k) ) / 2
//   u_w = ( u(i-1,j,k) + u(i-1,j+1,k) ) / 2
//


namespace Tuna {

template<class Tprec, int Dim>
inline bool CDS_ZCoDi<Tprec, Dim>::calcCoefficients3D () {
    prec_t dyz = dy * dz, dyz_dx = Gamma * dyz / dx;
    prec_t dxz = dx * dz, dxz_dy = Gamma * dxz / dy;
    prec_t dxy = dx * dy, dxy_dz = Gamma * dxy / dz;
    prec_t dxyz_dt = dx * dy * dz / dt;
    prec_t ce, cw;
    prec_t cn, cs;
    prec_t cf, cb;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int k = bk; k <= ek; ++k)
       for (int i =  bi; i <= ei; ++i)
	   for (int j = bj; j <= ej; ++j)
	   {
	       ce = ( u(i,j,k) + u(i,j+1,k) ) * 0.5 * dyz;
	       cw = ( u(i-1,j,k) + u(i-1,j+1,k) ) * 0.5 * dyz;
	       cn = ( v(i,j,k) + v(i+1,j,k) ) * 0.5 * dxz;
	       cs = ( v(i,j-1,k) + v(i+1,j-1,k) ) * 0.5 * dxz;
	       cf = ( w(i,j,k) + w(i,j,k+1) ) * 0.5 * dxy;
	       cb = ( w(i,j,k) + w(i,j,k-1) ) * 0.5 * dxy;

	       aE (i,j,k) = (dyz_dx - ce * 0.5);
	       aW (i,j,k) = (dyz_dx + cw * 0.5);
	       aN (i,j,k) = (dxz_dy - cn * 0.5);
	       aS (i,j,k) = (dxz_dy + cs * 0.5);
	       aF (i,j,k) = (dxy_dz - cf * 0.5);
	       aB (i,j,k) = (dxy_dz + cb * 0.5);
	       aP (i,j,k) = aE (i,j,k) + aW (i,j,k) + aN (i,j,k) + aS (i,j,k)
		   + aF (i,j,k) + aB (i,j,k) + dxyz_dt;
//		+ (ce - cw) + (cn - cs) + (cf - cb);
// Term (ce - cw) is part of discretizated continuity equation, and
// must be equal to zero when that equation is valid, so I can avoid
// this term for efficiency.
	       sp(i,j,k) += w(i,j,k) * dxyz_dt - 
		 ( p(i,j,k+1)- p(i,j,k) ) * dxy; 
	   }
   calc_dw_3D();
   applyBoundaryConditions3D();
   return 0;
}


} // Tuna namespace













